PairDirection Subroutine

SUBROUTINE PairDirection &
!
(stations)


IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: stations

!Local declarations
  INTEGER (KIND = short) :: i, j
  
!---------------------------end of declarations--------------------------------

  
!allocate variables
ALLOCATE ( dir_pairs (stations % countObs,stations % countObs) )

! Compute direction angle between stations 
dir_pairs = 0.
DO i = 1, stations % countObs
    DO j = 1, stations % countObs
        
        IF ( i == j) CYCLE !skip same point pairs, distance is zero
        
        IF ( j < i) CYCLE !this pair is already in the vector; pair(4,2) is the same as pair(2,4)
        
        point1 % northing = stations % obs (i) % xyz % northing  
        point1 % easting = stations % obs (i) % xyz % easting 
          
        point2 % northing = stations % obs (j) % xyz % northing  
        point2 % easting = stations % obs (j) % xyz % easting
                
        dir_pairs(i,j) = DirectionAngle (point1,point2)   
        
        !direction angle is measured from North clockwise.
        !change convention to angle measured from x- axis (east) anticlockwise.
        IF ( dir_pairs (i,j) <= pi / 2. ) THEN
            dir_pairs (i,j) = pi /2. - dir_pairs (i,j)
        ELSE
            !dir_pairs (i,j) = 3./2.*pi - dir_pairs (i,j)
            dir_pairs (i,j) =   pi + dir_pairs (i,j)
        END IF
       
    END DO
END DO 

RETURN
END SUBROUTINE PairDirection